Mapping

Download mapping.Rmd

Download mapping.R

Common solutions for creating maps usually involves GIS software, such as ArcGIS, QGIS, eSpatial, etc., which allow to visually prepare a map, in the same approach as one would prepare a poster or a document layout. On the other hand, R has also advanced spatial capabilities and can be used to draw maps. As we mentioned earlier for the graphics, R allows you to create map highly customization which makes it a solution increasingly attractive.

Several packahges will be used today. It is generally a good practice to group the package you will use at the beginning of your script:

library(maptools)
library(rgdal)
library (sp)
library(raster)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
library(rgbif)
library(mapr)
library(marmap)
library(leaflet)

#Sys.setlocale("LC_TIME", "English")
#Sys.setlocale(category = "LC_ALL", "Chinese (Traditional)_Taiwan.950")

Simple map

Getting a simple map is pretty straightforward in R using the packagemaptools

data(wrld_simpl) 
# NOT RUN: ? wrld_simpl
plot(wrld_simpl)

Now we know the options to modify a plot, therefore you can customize this map with some elements we mentioned earlier:

plot(wrld_simpl,xlim=c(100,130),ylim=c(-20,50),col='olivedrab3',bg='lightblue')

Try to zoom on Taiwan: you will quickly realize that we have a pretty bad resolution a finer scale, see ?wrld_simpl

Shape files

If you like running, hiking, walking, or being lazy in a park - many devices now can export gpx format. You an read the file directly as lines (i.e, tracks), points (i.e, track_points), and a few other formats you can find in the help for readOGR from the famous package rgdal: bindings for the ‘Geospatial’ Data Abstraction Library. Note that the same could be accomplished using the sf package and st_readfunction.

par(mfrow=c(1,2))
run1 <- readOGR(dsn="Data/run.gpx",layer="tracks")
OGR data source with driver: GPX 
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data\run.gpx", layer: "tracks"
with 1 features
It has 12 fields
plot(run1, main='Line') # my running activity
run2 <- readOGR(dsn="Data/run.gpx",layer="track_points")
OGR data source with driver: GPX 
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data\run.gpx", layer: "track_points"
with 931 features
It has 26 fields
plot(run2, main='Points')
dev.off()
RStudioGD 
        2

The reverse can be accomplish with the function writeOGR to export objects as SpatialPointsDataFrame, SpatialLinesDataFrame, or SpatialPolygonsDataFrame (as defined in the sp package). This will write an ArcGIS compatible shapefile, but many different formats are available by specifying the correct driver.Now that we know how to get a basic map in R, let’s look at how we can export and import data.

writeOGR(wrld_simpl, dsn="Data", layer = "world_test", driver = "ESRI Shapefile", overwrite_layer = TRUE)

Now we could open world_test.shp in ArcGIS (or others), but we can also import shapefile (.shp) back into R, let’s use that same file

world_shp <- readOGR(dsn = "Data",layer = "world_test")
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data", layer: "world_test"
with 246 features
It has 11 fields
plot(world_shp)

Spatial data types

Creating spatial data from scratch in R seems a little bit confusing at first, but once you got the logic behind, it get easier. As we learned before for graphics you can add vector based object (points, Lines, polygones) to your map. A transformation of the numerical x and y values applied in this case to recognize those values as spatial coordinate.

plot(wrld_simpl,xlim=c(115,128) ,ylim=c(19.5,27.5),col='#D2B48C',bg='lightblue') # TW map
coords <- matrix(c(121.537290,120.265541, 25.021335, 22.626524),ncol=2) # NTU and SYS univ. 
coords <- coordinates(coords) # assign values as spatial coordinates
spoints <- SpatialPoints(coords) # create SpatialPoints
df <- data.frame(location=c("NTU","SYS")) # create a dataframe
spointsdf <- SpatialPointsDataFrame(spoints,df) # create a SpatialPointsDataFrame
plot(spointsdf,add=T,col=c('black','black'),pch=19,cex=2.2) # plot it on our map
text(121,24, 'TAIWAN', cex=1)

plot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
l <- Line(coords)
ls <- Lines(list(l),ID="1")
sls <- SpatialLines(list(ls))
df <- data.frame(province="Saskatchewan")
sldf <- SpatialLinesDataFrame(sls,df)
plot(sldf,add=T,col='#3d2402', cex=2)
text(-114, 55, 'Saskatchewan', srt=90, cex=0.7)
text(-114, 63, 'CANADA', cex=1)

plot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
p <- Polygon(coords)
ps <- Polygons(list(p),ID="1")
sps <- SpatialPolygons(list(ps))
df <- data.frame(province="Saskatchewan")
spdf <- SpatialPolygonsDataFrame(sps,df)
plot(spdf,add=T,col='#45220d') 
text(-114, 55, 'Saskatchewan', srt=90, cex=0.7)
text(-114, 63, 'CANADA', cex=1)
text(-103, 46, 'UNITED STATES', cex=1)
text(-40, 78, 'GREENLAND', cex=1)
text(-35, 55, 'Atlantic Ocean', cex=1, col='#071238')

raster fct.

Spatial data: download, unzip, and import spatial shape (don’t forget to set up your working directory). Those spatial data are high resolution spatial data at the country level.

Using the rasterpackage and the function getData, you can easily download directly polygones (vectorized) shape from GADM (Global Administrative Areas, or other sources). You will simply need the country code such as TWN for Taiwan, JPN for Japan, etc. In addition, the argument level indicates the level you are trying to get. For provinces (first level subdivision)) level=1, for the overall country level level=0. Check ?getDta.

TWN1 <- getData('GADM', country="TWN", level=0) # data Taiwan
JPN <- getData('GADM', country="JPN", level=0) # data Japan
class(TWN1) # those datasets are SpatialPolygonsDataFrame
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
par(mfrow = c(1, 2))
plot(TWN1,axes=T,bg=colors()[431],col='grey')
plot(JPN,axes=T,bg=colors()[431],col='grey')

Don’t forget to close your graphical window:

dev.off()
RStudioGD 
        2

Simply set up xlim and ylim

plot (TWN1, axes=T, xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431],col='grey') 

TWN2 <- getData('GADM', country="TWN", level=1)
TWN2$NAME_1
[1] "Fujian"     "Kaohsiung"  "New Taipei" "Taichung"   "Tainan"    
[6] "Taipei"     "Taiwan"

You can see the list of counties is not complete. It is often the case when using GADM that some ‘details’ will be missing. Let’s visualize one county for which we have the data shape.

plot(TWN1,col="grey",xlim=c(119,122.5), ylim=c(21.5,25.5), bg=colors()[431], axes=T)
KAO <- TWN2[TWN2$NAME_1=="Kaohsiung",]
plot(KAO,col="grey 33",add=TRUE)

Note the use of add=TRUE in the latest plot function.

# base map
plot(TWN1,col="grey",xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431], axes=T)
# adding  spatial polygones 
TAI <- TWN2[TWN2$NAME_1=="Taipei" | TWN2$NAME_1=="New Taipei" ,]
plot(TAI,col="black",add=TRUE)
# adding spatial points 
coords <- matrix(cbind(lon=c(121.2,121.55,121.8),lat=c(25,25.19,24.5)),ncol=2)
coords <- coordinates(coords)
spoints <- SpatialPoints(coords)
df <- data.frame(location=c("City 1","City 2","City 3"),pop=c(138644,390095,34562))
spointsdf <- SpatialPointsDataFrame(spoints,df)
scalefactor <- sqrt(spointsdf$pop)/sqrt(max(spointsdf$pop))
plot(spointsdf,add=TRUE,col='white',pch=1,cex=scalefactor*3,lwd=2) 
# adding a location of NTU (not spatial point here)
points(121.537290,25.021335, type="p", pch=18, col='white', cex=1.5)
# adding text
text(121.53,24.921335,"NTU", col='white', font=2)
# adding scale
maps::map.scale(x=120, y=25.4)
# adding north arrow
GISTools::north.arrow(xb=120.3,yb=24.7, len=0.06, lab='N')

Note the use of the synthax maps::map.scale. It means you wanna use the function map.scale from the package maps. It avoid passing by library(maps), thus avoiding possible conflict among functions.

Each country will usually have open data platform with country shape data. Those data are usually more accurate and should be preferred. In Taiwan, shapefile of counties are available for download here. You could get even get a finer resolution at the township level if you wish (large file).

First prepare your system to receive traditional Chinese characters:

Sys.setlocale(category = "LC_ALL", "Chinese (Traditional)_Taiwan.950")
[1] "LC_COLLATE=Chinese (Traditional)_Taiwan.950;LC_CTYPE=Chinese (Traditional)_Taiwan.950;LC_MONETARY=Chinese (Traditional)_Taiwan.950;LC_NUMERIC=C;LC_TIME=Chinese (Traditional)_Taiwan.950"
sessionInfo()
R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Traditional)_Taiwan.950 
[2] LC_CTYPE=Chinese (Traditional)_Taiwan.950   
[3] LC_MONETARY=Chinese (Traditional)_Taiwan.950
[4] LC_NUMERIC=C                                
[5] LC_TIME=Chinese (Traditional)_Taiwan.950    
system code page: 65001

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] psych_2.2.5             MVPARTwrap_0.1-9.2     
 [3] mvpart_1.6-2            mvabund_4.2.1          
 [5] devtools_2.4.4          usethis_2.1.6          
 [7] caret_6.0-93            randomForest_4.7-1.1   
 [9] rattle_5.5.1            bitops_1.0-7           
[11] tibble_3.1.7            rpart_4.1.16           
[13] tree_1.0-42             e1071_1.7-11           
[15] gridExtra_2.3           fpc_2.2-9              
[17] factoextra_1.0.7        qgraph_1.9.2           
[19] gclus_1.3.2             cluster_2.1.3          
[21] ade4_1.7-19             ape_5.6-2              
[23] vegan_2.6-2             permute_0.9-7          
[25] marmap_1.0.6            mapr_0.5.2             
[27] rgbif_3.7.2             ggspatial_1.1.6        
[29] rnaturalearthdata_0.1.0 rnaturalearth_0.1.0    
[31] sf_1.0-8                raster_3.5-29          
[33] rgdal_1.5-32            maptools_1.1-4         
[35] sp_1.5-0                gganimate_1.0.7        
[37] animation_2.7           nlme_3.1-157           
[39] lmerTest_3.1-3          lme4_1.1-30            
[41] readr_2.1.2             yarrr_0.1.5            
[43] circlize_0.4.15         BayesFactor_0.9.12-4.4 
[45] Matrix_1.4-1            coda_0.19-4            
[47] jpeg_0.1-9              interactions_1.1.5     
[49] car_3.1-0               carData_3.0-5          
[51] MASS_7.3-57             corrplot_0.92          
[53] Hmisc_4.7-1             Formula_1.2-4          
[55] survival_3.3-1          leaflet_2.1.1          
[57] forcats_0.5.1           hexbin_1.28.2          
[59] stringr_1.4.0           tidyr_1.2.0            
[61] ggplot2_3.3.6           lattice_0.20-45        
[63] readxl_1.4.0            dplyr_1.0.9            
[65] scales_1.2.0           

loaded via a namespace (and not attached):
  [1] corpcor_1.6.10          class_7.3-20           
  [3] ps_1.7.1                foreach_1.5.2          
  [5] rprojroot_2.0.3         crayon_1.5.1           
  [7] fBasics_4021.92         backports_1.4.1        
  [9] rlang_1.0.3             nloptr_2.0.3           
 [11] callr_3.7.0             bit64_4.0.5            
 [13] glue_1.6.2              parallel_4.2.1         
 [15] processx_3.6.1          classInt_0.4-7         
 [17] tidyselect_1.1.2        ggpubr_0.4.0           
 [19] xtable_1.8-4            MatrixModels_0.5-0     
 [21] magrittr_2.0.3          evaluate_0.16          
 [23] ncdf4_1.19              cli_3.3.0              
 [25] rstudioapi_0.13         miniUI_0.1.1.1         
 [27] whisker_0.4             bslib_0.3.1            
 [29] jtools_2.2.0            maps_3.4.0             
 [31] GISTools_0.7-4          shiny_1.7.2            
 [33] xfun_0.31               pkgbuild_1.3.1         
 [35] ggrepel_0.9.1           listenv_0.8.0          
 [37] installr_0.23.2         png_0.1-7              
 [39] future_1.27.0           ipred_0.9-13           
 [41] withr_2.5.0             tweedie_2.3.5          
 [43] plyr_1.8.7              cellranger_1.1.0       
 [45] hardhat_1.2.0           pROC_1.18.0            
 [47] pillar_1.8.1            GlobalOptions_0.1.2    
 [49] cachem_1.0.6            fs_1.5.2               
 [51] flexmix_2.3-18          kernlab_0.9-31         
 [53] vctrs_0.4.1             pbivnorm_0.6.0         
 [55] ellipsis_0.3.2          oai_0.3.2              
 [57] generics_0.1.3          lava_1.6.10            
 [59] urltools_1.7.3          tools_4.2.1            
 [61] foreign_0.8-82          munsell_0.5.0          
 [63] tweenr_2.0.0            proxy_0.4-27           
 [65] fastmap_1.1.0           compiler_4.2.1         
 [67] pkgload_1.3.0           abind_1.4-5            
 [69] httpuv_1.6.5            sessioninfo_1.2.2      
 [71] rgeos_0.5-9             prodlim_2019.11.13     
 [73] rpart.plot_3.1.1        deldir_1.0-6           
 [75] utf8_1.2.2              later_1.3.0            
 [77] recipes_1.0.1           wk_0.6.0               
 [79] jsonlite_1.8.0          pbapply_1.5-0          
 [81] lazyeval_0.2.2          promises_1.2.0.1       
 [83] latticeExtra_0.6-30     checkmate_2.1.0        
 [85] rmarkdown_2.14          textshaping_0.3.6      
 [87] statmod_1.4.37          Rtsne_0.16             
 [89] pander_0.6.5            igraph_1.3.4           
 [91] numDeriv_2016.8-1.1     yaml_2.3.5             
 [93] systemfonts_1.0.4       prabclus_2.3-2         
 [95] conditionz_0.1.0        htmltools_0.5.2        
 [97] memoise_2.0.1           lavaan_0.6-12          
 [99] profvis_0.3.7           modeltools_0.2-23      
[101] viridisLite_0.4.0       digest_0.6.29          
[103] assertthat_0.2.1        mime_0.12              
[105] adehabitatMA_0.3.14     units_0.8-0            
[107] RSQLite_2.2.16          future.apply_1.9.0     
[109] remotes_2.4.2           data.table_1.14.2      
[111] urlchecker_1.0.1        blob_1.2.3             
[113] ragg_1.2.2              ggsci_2.9              
[115] splines_4.2.1           labeling_0.4.2         
[117] broom_1.0.0             hms_1.1.1              
[119] colorspace_2.0-3        base64enc_0.1-3        
[121] mnormt_2.1.0            shape_1.4.6            
[123] nnet_7.3-17             sass_0.4.1             
[125] Rcpp_1.0.9              mclust_5.4.10          
[127] bookdown_0.28           mvtnorm_1.1-3          
[129] fansi_1.0.3             tzdb_0.3.0             
[131] parallelly_1.32.1       ModelMetrics_1.2.2.2   
[133] R6_2.5.1                grid_4.2.1             
[135] crul_1.2.0              lifecycle_1.0.1        
[137] formatR_1.12            zip_2.2.0              
[139] ggsignif_0.6.3          curl_4.3.2             
[141] minqa_1.2.4             jquerylib_0.1.4        
[143] robustbase_0.95-0       distill_1.4            
[145] glasso_1.11             RColorBrewer_1.1-3     
[147] iterators_1.0.14        gower_1.0.0            
[149] htmlwidgets_1.5.4       triebeard_0.3.0        
[151] purrr_0.3.4             crosstalk_1.2.0        
[153] terra_1.6-7             labdsv_2.0-1           
[155] timeSeries_4021.104     mgcv_1.8-40            
[157] globals_0.16.0          htmlTable_2.4.1        
[159] codetools_0.2-18        lubridate_1.8.0        
[161] gtools_3.9.3            leaflet.providers_1.9.0
[163] prettyunits_1.1.1       gtable_0.3.0           
[165] DBI_1.1.3               stats4_4.2.1           
[167] httr_1.4.3              highr_0.9              
[169] KernSmooth_2.23-20      stringi_1.7.6          
[171] progress_1.2.2          spatial_7.3-15         
[173] reshape2_1.4.4          farver_2.1.1           
[175] uuid_1.1-0              diptest_0.76-0         
[177] fdrtool_1.2.17          magick_2.7.3           
[179] timeDate_4021.104       xml2_1.3.3             
[181] boot_1.3-28             s2_1.1.0               
[183] interp_1.1-3            DEoptimR_1.0-11        
[185] bit_4.0.4               pkgconfig_2.0.3        
[187] rstatix_0.7.0           knitr_1.39             
[189] downlit_0.4.2           httpcode_0.3.0

You can go to the website and downaload data. Alternatively, using the code source of the page your can target directly the the url of teh dataset you want to download. This is a zip file that you want to temporary store and extract in a specific directory

url <- 'https://data.moi.gov.tw/MoiOD/System/DownloadFile.aspx?DATA=72874C55-884D-4CEA-B7D6-F60B0BE85AB0'
path1 <- tempfile(fileext = ".zip")
if (file.exists(path1))  'file alredy exists' else download.file(url, path1, mode="wb")
zip::unzip(zipfile = path1,exdir = 'Data')

Again, the function readOGR reads the .shp file. Other arguments are for the file encoding and to deal with the Chinese characters because my system is in English.

taiwan <- readOGR('Data/COUNTY_MOI_1090820.shp', use_iconv=TRUE, encoding='UTF-8')
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data\COUNTY_MOI_1090820.shp", layer: "COUNTY_MOI_1090820"
with 22 features
It has 4 fields

Make the plot:

plot(taiwan,axes=T,bg=colors()[431],col='grey')

This is the official map of TAIWAN, including all remote territories such as Taiping and Diaoyutai islands. County names are here accessible using:

taiwan$COUNTYNAME
 [1] "連江縣" "宜蘭縣" "彰化縣" "南投縣" "雲林縣" "基隆市" "臺北市"
 [8] "新北市" "臺中市" "臺南市" "桃園市" "苗栗縣" "嘉義市" "嘉義縣"
[15] "金門縣" "高雄市" "臺東縣" "花蓮縣" "澎湖縣" "新竹市" "新竹縣"
[22] "屏東縣"
taiwan$COUNTYENG 
 [1] "Lienchiang County" "Yilan County"      "Changhua County"  
 [4] "Nantou County"     "Yunlin County"     "Keelung City"     
 [7] "Taipei City"       "New Taipei City"   "Taichung City"    
[10] "Tainan City"       "Taoyuan City"      "Miaoli County"    
[13] "Chiayi City"       "Chiayi County"     "Kinmen County"    
[16] "Kaohsiung City"    "Taitung County"    "Hualien County"   
[19] "Penghu County"     "Hsinchu City"      "Hsinchu County"   
[22] "Pingtung County"

ggplot2maps

As we mentioned earlier, the package ggplot2 implements the grammar of graphics in R. While ggplot2 is becoming the standard for R graphs, it does not handle spatial data specifically. The current state-of-the-art of spatial objects in R relies on Spatial classes defined in the package sp, but the new package sf has implemented the “simple feature” standard, and is steadily taking over sp. Recently, the package ggplot2 has allowed the use of simple features from the package sf as layers in a graph (since the version 3.0.0 of ggplot2). The combination of ggplot2 and sf therefore enables to programmatically create maps, using the grammar of graphics, just as informative or visually appealing as traditional GIS software.

Note that the conversion from sp to sf is acheivable using the function st_as_sf from the package sf.

# not run
# wrld_simpl <- st_as_sf(wrld_simpl)

Getting started: theme and datasets

After loading the basic packages necessary for all maps, i.e. ggplot2 and sf. We can change the theme of our plot. The classic dark-on-light theme for ggplot2 (theme_bw) is appropriate for maps (default: theme_grey()):

theme_set(theme_bw())

The package rnaturalearth provides a map of countries of the entire world. Use ne_countries to pull country data and choose the scale (rnaturalearthhires is necessary for scale = "large"). The function can return sp classes (default) or directly sf classes, as defined in the argument returnclass:

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
[1] "sf"         "data.frame"

In addition to country shape polygones, this dataset contains information on the population in every country.

ggplot map

First, let us start with creating a base map of the world using ggplot2. This base map will then be extended with different map elements, as well as zoomed in to an area of interest. We can check that the world map was properly retrieved and converted into an sf object, and plot it with ggplot2. This time we will indicate that the geometry of our plot is constrained by a sf object:

ggplot(data = world) +
   geom_sf()

In your map of the world, the plot panel is expanded beyond the size of the earth (you can see that the graticule lines end before the edge of the plot panel), and hence no axis ticks are drawn. One way to solve the issue is to turn off the expansion while defining the coordinates.

ggplot(data = world) +
   geom_sf() +
   coord_sf(expand = FALSE)

As before, the layers are added one at a time in a ggplot call, so the order of each layer is very important. All data will have to be in an sf format to be used by ggplot2; data in other formats (e.g. classes from sp) will be manually converted to sf classes if necessary.

A title and a subtitle can be added to the map using the function ggtitle, passing any valid character string (e.g. with quotation marks) as arguments. Axis names are absent by default on a map, but can be changed to something more suitable (e.g. “Longitude” and “Latitude”), depending on the map:

ggplot(data = world) +
    geom_sf() +
    coord_sf(expand = FALSE) +
    xlab("Longitude") + ylab("Latitude") +
    ggtitle("World map", subtitle = paste0("(", length(unique(world$name)), " countries)"))

In many ways, sf geometries are no different than regular geometries, and can be displayed with the same level of control on their attributes. Here is an example with the polygons of the countries filled with a green color (argument fill), using black for the outline of the countries (argument color):

ggplot(data = world) + 
    geom_sf(color = "black", fill = "lightgreen") +
    coord_sf(expand = FALSE)

The package ggplot2 allows the use of more complex color schemes, such as a gradient on one variable of the data. Here is another example that shows the population of each country. In this example, we use the “viridis” colorblind-friendly palette for the color gradient (with option = "magma" for the magma variant), using the square root of the population (which is stored in the variable POP_EST of the world object):

ggplot(data = world) +
    geom_sf(aes(fill = pop_est)) +
    coord_sf(expand = FALSE) +
    scale_fill_viridis_c(option = "plasma", trans = "sqrt") # sqrt transform

The function coord_sf allows to deal with the coordinate system, which includes both projection and extent of the map. By default, the map will use the coordinate system of the first layer that defines one (i.e. scanned in the order provided), or if none, fall back on WGS84 (latitude/longitude, the reference system used in GPS). Using the argument crs, it is possible to override this setting, and project on the fly to any projection. This can be achieved using any valid PROJ4 string (here, the European-centric ETRS89 Lambert Azimuthal Equal-Area projection):

ggplot(data = world) +
    geom_sf() +
    scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
    coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

Spatial Reference System Identifier (SRID) or an European Petroleum Survey Group (EPSG) code are available for the projection of interest, they can be used directly instead of the full PROJ4 string. The two following calls are equivalent for the ETRS89 Lambert Azimuthal Equal-Area projection, which is EPSG code 3035:

ggplot(data = world) +
    geom_sf() +
    coord_sf(crs = "+init=epsg:3035")

# OR

ggplot(data = world) +
    geom_sf() +
    coord_sf(crs = st_crs(3035))

The extent of the map can also be set in coord_sf, in practice allowing to “zoom” in the area of interest, provided by limits on the x-axis (xlim), and on the y-axis (ylim). Note that the limits are automatically expanded by a fraction to ensure that data and axes don’t overlap; it can also be turned off to exactly match the limits provided with expand = FALSE:

ggplot(data = world) +
    geom_sf(aes(fill = pop_est)) +
    coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE) +
    scale_fill_viridis_c(option = "plasma") # linear scale

ggspatial elements

Several packages are available to create a scale bar on a map (e.g. prettymapr, vcd, ggsn, or legendMap). Here the package ggspatial provides easy-to-use functions. scale_bar that allows to add simultaneously the north symbol and a scale bar into the ggplot map. Five arguments need to be set manually: lon, lat, distance_lon, distance_lat, and distance_legend. The location of the scale bar has to be specified in longitude/latitude in the lon and lat arguments. The shaded distance inside the scale bar is controlled by the distance_lon argument while its width is determined by distance_lat. Additionally, it is possible to change the font size for the legend of the scale bar (argument legend_size, which defaults to 3). The North arrow behind the “N” north symbol can also be adjusted for its length (arrow_length), its distance to the scale (arrow_distance), or the size the N north symbol itself (arrow_north_size, which defaults to 6). Note that all distances (distance_lon, distance_lat, distance_legend, arrow_length, arrow_distance) are set to "km" by default in distance_unit; they can also be set to nautical miles with “nm”, or miles with “mi”.

ggplot(data = world) +
    geom_sf(aes(fill = pop_est)) +
    coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE) +
    scale_fill_viridis_c(option = "plasma") +
    annotation_scale(location = "br", width_hint = 0.3) +
    annotation_north_arrow(location = "br", which_north = "true", 
        pad_x = unit(0.5, "cm"), pad_y = unit(1, "cm"),
        style = north_arrow_fancy_orienteering)

Note: if you plot a larger area, you may get a warning on the inaccuracy of the scale bar.

The world data set already contains country names and the coordinates of the centroid of each country (among more information). We can use this information to plot country names, using world as a regular data.frame in ggplot2. The function geom_text can be used to add a layer of text to a map using geographic coordinates. The function requires the data needed to enter the country names, which is the same data as the world map. Again, we have a very flexible control to adjust the text at will on many aspects:

Additionally, the annotate function can be used to add a single character string at a specific location, as demonstrated here to add the “Gulf of Mexico”Pacific Ocean” and Ryukyu Archipelago

sf::sf_use_s2(FALSE) # FOR ERROR turn off the s2 processing
world_points<- st_centroid(world)
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))

ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
geom_text(data= world_points,aes(x=X, y=Y, label=name),
    color = "grey", fontface = "bold", check_overlap = FALSE) +
    scale_fill_viridis_c(option = "plasma") +
    annotate(geom = "text", x = 124, y = 21, label = "Pacific Ocean", fontface = "italic", color = "#0b3c8a", size = 5) +
   annotate(geom = "text", x = 124.2, y = 24, label = "Ryukyu archipelago", fontface = "italic", color = "#d41919", size = 3) + 
   coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE)

To goo further: Finalization

Now to make the final touches, the theme of the map can be edited to make it more appealing. We suggested the use of theme_bw for a standard theme, but there are many other themes that can be selected from (see for instance ?ggtheme in ggplot2, or the package ggthemes which provide several useful themes). Moreover, specific theme elements can be tweaked to get to the final outcome:

ggplot(data = world) +
   geom_sf(fill= 'antiquewhite') + 
   geom_text(data= world_points,aes(x=X, y=Y, label=name), color = 'darkblue', fontface = 'bold', check_overlap = FALSE) + 
   annotate(geom = 'text', x = -90, y = 26, label = 'Gulf of Mexico', fontface = 'italic', color = 'grey22', size = 6) +
   annotation_scale(location = 'bl', width_hint = 0.5) + 
   annotation_north_arrow(location = 'bl', which_north = 'true', pad_x = unit(0.75, 'in'), pad_y = unit(0.5, 'in'), style = north_arrow_fancy_orienteering) + 
   coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE) + 
   xlab('Longitude') + ylab('Latitude') + 
   ggtitle('Map of the Gulf of Mexico and the Caribbean Sea') + 
   theme(panel.grid.major = element_line(color = gray(.5), linetype = 'dashed', size = 0.5), panel.background = element_rect(fill = 'aliceblue'))

Exporting map

The final map is now ready, it is very easy to save it using ggsave. This function allows a graphic (typically the last plot displayed) to be saved in a variety of formats, including the most common PNG (raster bitmap) and PDF (vector graphics), with control over the size and resolution of the outcome. For instance here, we save a PDF version of the map, which keeps the best quality, and a PNG version of it for web purposes:

ggsave("Figures/Datamap.pdf")
ggsave("Figures/map_web.png", width = 6, height = 6, dpi = "screen")

Special maps

Species distribution rgbif

Import Global Biodiversity Information Facility (GBIF) data using rgbif package and crete distribution map using the mapr package.

gbif.res <- occ_search(scientificName = "Urocissa caerulea", limit = 1200)

map_ggplot(gbif.res) +
  coord_sf(xlim = c(120, 123), ylim = c(21, 26), expand = FALSE)

Note the slight difference in overlap of the distribution data on this shape file.

Have a look here, to combine climate and species data.

Bathymetric map with marmap

marmap can query and plot NOAA’s bathymetry database

# query
TW.bathy <- getNOAA.bathy(lon1=118,lon2=124, lat1=21,lat2=26,resolution=1) # don't put too wide / resolution: 1 
# define palette
blues <- colorRampPalette(c("darkblue", "cyan"))
greys <- colorRampPalette(c(grey(0.4),grey(0.99)))
# make the plot
plot.bathy(TW.bathy,
     image=T,
     deepest.isobath=c(-6000,-120,-30,0),
     shallowest.isobath=c(-1000,-60,0,0),
     step=c(2000,60,30,0),
     lwd=c(0.3,1,1,2),
     lty=c(1,1,1,1),
     col=c("grey","black","black","black"), 
     drawlabels=c(T,T,T,F),
     bpal = list(c(0,max(TW.bathy),greys(100)),c(min(TW.bathy),0,blues(100))),
     land=T, xaxs="i"
     )

Profiles can be extract using get.transect:

tw.profile <-get.transect(TW.bathy,x1=119.5,y1=23.75, x2=122,y2=23.75, dis=TRUE)
plotProfile(tw.profile) 
#### Not Run: extract a profile Manually
#### manual.profile<-get.transect (TW.bathy, loc=T,dist=T) 
#### plotProfile(manual.profile)

Interactive maps

Leaflet is one of the most popular open-source JavaScript libraries for interactive maps. You can loose hours and its functionalities are amazing. The R package leaflet makes it easy to integrate and control Leaflet maps in R.

FRE <- paste(sep = "<br/>",
  "<b><a href='https://www.dipintothereef.com/'>FRELAb TAIWAN</a></b>",
  "Functional Reef Ecology Lab",
  "Institute of Oceanography, NTU")


leaflet(taiwan) %>%
  addPolygons(weight=0.5) %>%
  addTiles(group="Kort") %>%
  addPopups(121.53725, 25.021252, FRE, options = popupOptions(closeButton = FALSE))


⚠ Practice 3: Create and customize an interactive map of your choice after exploring the functionality of the leaflet package. Push both your .Rmd and .html files into a public repository available from your Github account.You will share with me be email [] the address (URL) of this repository (such as https://github.com/vianneydenis/OCEAN5098.git) before next Monday in order for me to check your work. The title of your email should be `Practice 3 (your name: your student no.). BE CREATIVE ;)

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".